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Abstract 

We illustrate the implementation of a method based on the use of recursion 
relations in (Bjorken) x— space for the solution of the evolution equations of QCD 
for all the leading twist distributions. The algorithm has the advantage of being 
very fast. The implementation that we release is written in C and is performed to 
next-to-leading order in a s . 
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Program summary 

• Title of program: evolution, c 

• Computer: Athlon 1800 plus 

• Operating system under which the program has been tested: Linux 

• Programming language used: C 

• Peripherals used: Laser printer 

• No. of bytes in distributed program: 20771 

• Keywords: Polarized and unpolarized parton distribution, numerical solutions of 
renormalization group equation in Quantum Chromodynamics. 

• Nature of physical problem: The program provided here solves the DGLAP evolution 
equations to next-to-leading order a s , for unpolarized, longitudinally polarized and 
transversely polarized parton distributions. 

• Method of solution: We use a recursive method based on an expansion of the solution 
in powers of log (a s (Q) / a s (Q )) 

• Typical running time: About 1 minute and 30 seconds for the unpolarized and 
longitudinally polarized cases and 1 minute for the transversely polarized case. 

LONG WRITE UP 
1 Introduction 

Our understanding of the QCD dynamics has improved steadily along the years. In 
fact we can claim that precision studies of the QCD background in a variety of energy 
ranges, from intermediate to high energy - whenever perturbation theory and factorization 
theorems hold - are now possible and the level of accuracy reached in this area is due both 
to theoretical improvements and to the flexible numerical implementations of many of the 
algorithms developed at theory level. 

Beside their importance in the determination of the QCD background in the search of 
new physics, these studies convey an understanding of the structure of the nucleon from 
first principles, an effort which is engaging both theoretically and experimentally. 

It should be clear, however, that perturbative methods are still bound to approxima- 
tions, from the order of the perturbative expansion to the phenomenological description 
of the parton distribution functions, being them related to a parton model view of the 
QCD dynamics. 

Within the framework of the parton model, evolution equations of DGLAP-type - and 
the corresponding initial conditions on the parton distributions - are among the most 
important blocks which characterize the description of the quark-gluon interaction and, 
as such, deserve continuous attention. Other parts of this description require the compu- 
tation of hard scatterings with high accuracy and an understanding of the fragmentation 
region as well. The huge success of the parton model justifies all the effort. 

In this consolidated program, we believe that any attempt to study the renormalization 
group evolution describing the perturbative change of the distributions with energy, from 
a different - in this case, numerical - standpoint is of interest. 
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In this paper we illustrate an algorithm based on the use of recursion relations for 
the solution of evolution equations of DGLAP type. We are going to discuss some of the 
salient features of the method and illustrate its detailed implementation as a computer 
code up to next-to leading order (NLO) in a s , the strong coupling constant. 

In this context, it is worth to recall that the most common method implemented so 
far in the solution of the DGLAP equations is the one based on the Mellin moments, with 
all its good points and limitations. 

The reason for such limitations are quite obvious: while it is rather straightforward to 
solve the equations in the space of moments, their correct inversion is harder to perform, 
since this requires the computation of an integral in the complex plane and the search for 
an optimized path. 

In this respect, several alternative implementations of the NLO evolution are available 
from the previous literature, either based on the use of "brute force" algorithms [I] or on 
the Laguerre expansion 012], all with their positive features and their limitations. 

Here we document an implementation to NLO of a method based on an ansatz |3] 
which allows to rewrite the evolution equations as a set of recursion relations for some 
scale invariant functions, A n (x) and B n (x), which appear in the expansion. The advan- 
tage, compared to others, of using these recursion relations is that just few iterates of these 
are necessary in order to obtain a stable solution. The number of iterates is determined 
at run time. We also mention that our implementation can be extended to more com- 
plex cases, including the cases of nonforward parton distributions and of supersymmetry, 
which will be discussed elsewhere . Here we have implemented the recursion method in 
all the important cases of initial state scaling violations connected to the evolutions of 
both polarized and unpolarized parton densities, including the less known transverse spin 
distributions. 

Part of this document is supposed to illustrate the algorithm, having worked out in 
some detail the structure of the recursion relations rather explicitly, especially in the case 
of non-singlet evolutions, such as for transverse spin distributions. As we have mentioned, 
the same method has also been applied successfully to the case of supersymmetric QCD 
and will be illustrated in a companion paper. 

One of the advantages of the method is its analytical base, since the recursion relations 
can be written down in explicit form and at the same time one can perform a simple 
analytical matching between various regions in the evolutions, which is a good feature of 
x-spaced methods. While this last aspect is not relevant for ordinary QCD, it is relevant 
in other theories, such as for supersymmetric extensions of the parton model, once one 
assumes that, as the evolution scale Q raises, new anomalous dimensions are needed in 
order to describe the mixing between ordinary and supersymmetric partons. 



2 Definitions and Conventions 

In this section we briefly outline our definitions and conventions. 

We will be using the running of the coupling constant up to two-loop level 



a s (Q 2 ) 
where 
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and 

N 2 — 14 1 
Nc = 3, C F = ^- = -, Tf = T R n f = -n f , (3) 

where N c is the number of colors, n/ is the number of active flavors, which is fixed by the 
number of quarks with m q < Q. We have taken for the quark masses m c = 1.5 GeV, m h = 
4.5 GeV and m t = 175 GeV, these are necessary in order to identify the thresholds at 
which the number of flavours n./ is raised as we increase the final evolution scale. 

In our conventions Aqcd is denoted by A^^ and is given by 

A (3A5,6) = Q 24g ^ Q 20Q) Q 131) Q QeV _ ^ 

We also define the distribution of a given helicity (±), f ± (x,Q 2 ), which is the prob- 
ability of finding a parton of type / at a scale Q, where / = in a longitudinally 
polarized proton with the spin aligned (+) or anti-aligned (-) respect to the proton spin 
and carrying a fraction x of the proton's momentum. 

As usual, we introduce the longitudinally polarized parton distribution of the proton 

Af(x,Q 2 ) = f + (x,Q 2 )-f-(x,Q 2 ). (5) 

We also introduce another type of parton density, termed transverse spin distribution, 
which is defined as the probability of finding a parton of type / in a transversely polarized 
proton with its spin parallel (|) minus the probability of finding it antiparallel (J.) to the 
proton spin 

& T f(x,Q 2 ) = f(x,Q 2 ) - ft(x,Q 2 ). (6) 
Similarly, the unpolarized (spin averaged) parton distribution of the proton is given 

by 

f(x, Q 2 ) = Q 2 ) + f-(x, Q 2 ) = f( Xl Q 2 ) + f l (x, Q 2 ). (7) 

We also recall, if not obvious, that taking linear combinations of Equations (JZJ) and (jHJ, 
one recovers the parton distributions of a given helicity 

fy X ,Q>) = nX ' Q2)± 2 AnX ' Q2) . (8) 

In regard to the kernels, the notations P, AP, A T P, P +± , will be used to denote the 
Altarelli-Parisi kernels in the unpolarized, longitudinally polarized, transversely polarized, 
and the positive (negative) helicity cases respectively. 

Let us now turn to the evolution equations, starting from the unpolarized case. Defin- 
ing 



% (±) = ft ± ft, q {+) = E ^\ X* = # } " ~q {+ \ (9) 



i=l 



and introducing the convolution product 

f{x)®g{x)= [ 1( ^-f(-)g(y), (10) 



the evolution equations are 

d 



dlogQ 2 
d 



dlogQ 



qf\x,Q 2 ) = P NS -{x 1 a s {Q 2 ))®q ( f\x ) Q 2 ) ) (11) 
Xl (x,Q 2 ) = P NS+ (x,a s (Q 2 ))® X i(x,Q 2 ), (12) 
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for the non-singlet sector and 

p (n. ™ (r&\\ p ™ ('ns^ \ / „{+)(„ r&\ \ 

(13) 



d ( <?«(*, Q 2 ) \ ( P qq (x,a s (Q 2 )) P qg (x, a s (Q 2 )) \ ( q^{x,Q 2 ) 
dlogQ 2 I g(x,Q 2 ) ) I P gq (x,a s (Q 2 )) P gg (x } a s (Q 2 )) r 9 \ g(x,Q 2 ) 



for the singlet sector. Equations analogous to (CTTTni. with just a change of notation, 
are valid in the longitudinally polarized case and, due to the linearity of the evolution 
equations, also for the distributions in the helicity basis. In the transverse case instead, 
there is no coupling between gluons and quarks, so the singlet sector (fT3|) is missing. In 
this case we will solve just the nonsinglet equations 

d A T qt\x,Q 2 ) = A T P NS -(x,a s (Q 2 )) ® A T qt\x,Q 2 ), (14) 



dlogQ 2 



<J r Ar<?l +) (z,Q 2 ) = A T P NS+ (x,a s (Q 2 )) ® A T q^\x,Q 2 ). (15) 



dlogQ 2 

We also recall that the perturbative expansion, up to next-to-leading order, of the 
kernels is 

P(x, a.) = (g) P(°\x) + V>(*) + ■■■■ (16) 

Kernels of fixed helicity can also be introduced with P ++ (z) = (P(z) + AP(z))/2 and 
P + -(z) = (P(z) — AP(z))/2 denoting splitting functions of fixed helicity, which will be 
used below. 

The equations, in the helicity basis, are 

+Pl 9 + (-)®g4y) + Pl 9 -(-)®g-(y)), 
y y 

d -^f L = ^(P+-(-) ® g + (y) + ® 

at Itx y y 

+Pl 9 -{-)®g + {y) + Pl\{-)®g-{y)\ 
y y 

dt 2n y y 

+Pf+(-) (8 </ + (y) + Pf_(-) ® </-(*)), 

1/ y 
+PU-) ® 0+(v) + ® ^-(y)). (17) 

The non-singlet (valence) analogue of this equation is also easy to write down 

dj± ^ = ^(P ++ (~) ® ? +> v(y) + p+-(-) ® ?-,v(y)), 

dt In y y 

= ® ? + ,v(y) + ® ?-,v(y))- (18) 

dt In y y 

where the q^y = ?± — q± are the valence components of fixed helicity. The kernels in this 
basis are given by 

p(0) _ p(0) _ p (0) 

■ r NS±,++ ~ r qq,++ ~ r qq 
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= P w ,-+ = n f {x 


Pgq,++ 


= P qq,— = C F~ 


p(0) 


- P {0) — N ( 

- r gg,++ - iV c 1 


p(0) 
^99,+- 


= iV c ( 3x + - - 
V x 



I) 2 



- - 1 - x - x 2 ) + (3 Q 6(1 - x) 



(1 — x) + X 

3-x 2 ^. (19) 

Taking linear combinations of these equations (adding and subtracting) , one recovers 
the usual evolutions for unpolarized q(x) and longitudinally polarized Aq(x) distributions. 



3 The Ansatz and some Examples 

In order to solve the evolution equations directly in x-space, we assume solutions of the 
form 

f( x Q2) = y> An ^ \ Q „n a s{Q ) _i_ / Q 2n B n (x) a s (Q ) , v 

for each parton distribution /, where Q defines the initial evolution scale. The justifica- 
tion of this ansatz can be found, at least in the case of the photon structure function, in 
the original work of Rossi (4], and its connection to the ordinary solutions of the DGLAP 
equations is most easily worked out by taking moments of the scale invariant coefficient 
functions A n and B n and comparing them to the corresponding moments of the parton 
distributions, as we are going to illustrate in section 5. The link between Rossi's expansion 
and the solution of the evolution equations (which are ordinary differential equations) in 
the space of the moments up to NLO will be discussed in that section, from which it will 
be clear that Rossi's ansatz involves a resummation of the ordinary Mellin moments of 
the parton distributions. 

Setting Q = Q in J2Q|) we get 

f(x,Q 2 ) = A (x) + a s (Q 2 )B (x). (21) 

Inserting (|20|) in the evolution equations, we obtain the following recursion relations for 
the coefficients A n and B n 

A n+l {x) = -^P {0 \x) ® A n {x), (22) 
Po 

B n+1 (x) = -B n (x) - ^-A n+1 (x) - ^P {0) (x) (g) B n (x) - -^—P {l \x) ® A n (x), (23) 

obtained by equating left-hand sides and right-hand-side of the equation of the same 
logarithmic power in log n a s (Q 2 ) and a s log™ a s (Q 2 ). 

Any boundary condition satisfying (|2T|) can be chosen at the lowest scale Q and in 
our case we choose 

B (x) = 0, f(x,Q 2 )=A (x). (24) 

The actual implementation of the recursion relations is the main effort in the actual 
writing of the code. Obviously, this requires particular care in the handling of the singu- 
larities in x— space, being all the kernels defined as distributions. Since the distributions 



6 



are integrated, there are various ways to render the integrals finite, as discussed in the 
previous literature on the method [H] in the case of the photon structure function. In these 
previous studies the edge-point contributions - i.e. the terms which multiply 5(1 — x) in 
the kernels - are approximated using a sequence of functions converging to the 5 function 
in a distributional sense. 

This technique is not very efficient. We think that the best way to proceed is to 
actually perform the integrals explicitly in the recursion relations and let the subtracting 
terms appear under the same integral together with the bulk contributions (x < 1) (see 
also pj). This procedure is best exemplified by the integral relation 

dy fl , x f 1 dy yfjy) - xf(x) 

-f{x/y)= log(l - x)f(x) (25) 



y(l - y)+ Jx y y-x 

in which, on the right hand side, regularity of both the first and the second term is explicit. 
For instance, the evolution equations become (prior to separation between singlet and 
non-singlet sectors) in the unpolarized case 

d( li( x ) f dy yqi(y) - xqi(x) r 1 dy 3 



nn [ dy yqi(y - xq { (x f ^ dy 6 

d\o gm 2Cf J 7 — ^ — + 2CFlog{1 ~ xMx) ~ij {1+z) Qi{y) + 2 Cpq{x) 



C djL + (1 - zf) 9 (y) 
Jx y v 7 



+n f 

dg(x) n f l dy l + (l-z) 2 fi dy yf(y)-xf(x) 

Cf qi{y) + 2N c g(y) 



d\og(Q 2 ) Jx y z ' J x y y-x 

+2iV c log(l - x)g{x) + 2N C C (- - 2 + z(l - zj) g{y) + ^g(x) 

Jx y \z ) 2 

with z = x/y. Here q are fixed flavour distributions. 

4 An Example: The Evolution of the Transverse Spin 
Distributions 

Leading order (LO) and NLO recursion relations for the coefficients of the expansion can 
be worked out quite easily. We illustrate here in detail the implementation of a non-singlet 
evolutions, such as those involving transverse spin distributions. For the first recursion 
relation (f22|) in this case we have 

= -|"ArPW(x) ® At{x) = 



(26) 



C7 



A) 



C ^ yJM^M + At(x) log(1 _ x) 

Jx y y-x 



H£)tt d 7 A±M )M-B A±M - (27) 

As we move to NLO, it is convenient to summarize the structure of the transverse kernel 
ArP^ix) as 

A T P±' (1) (x) = K±{x)5{l -x) + K±(x)S 2 (x) + Kf{x) log(x) 

+Kf{x) log 2 (x) + Kf{x) \og{x) log(l -x) + Kf{x) 1 + K±{x) . (28) 

I 1 _ x )+ 
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Hence, for the (+) case we have 



ArP+^ix) ® A+(x) = K+A+(x) + £ ^ {K+(z)S 2 (z) + K 3 +(^) log(z) 



+ \og 2 (z)Kt(z) + log(z) log(l - z)K+{z)\ A+(y) + 

Kt \ f d ^y)-^ X ) + 1 ,1 dy 

[Jx y y-x J Jx y 

where z = x/y. For the (— ) case we get a similar expression. 
For the B^ +1 (x) we get (for the (+) case) 



2CV 



dyyA+(y) -xA+(x) 



x y 



y-x 



+ A+(rr) log(l - x) 



+ 



+ X 3 + (z) log(z) + log 2 (z)K 4 + (z) + log(z) log(l - z)Kt(z) 



+ 



i dyyA+(y) - xA+{x) 



x y y-x 
^dyyBt{y)-xBt{x) 



+ A+(x) log(l-x) 



+ K 7 f ~ A n(y) 

Jx y 



x y 



y-x 



+ B±(x)\og(l-x) 



+ 



where in the (+) case we have the expressions 



K+(x) = ^-C F (-2n f (3 + 4vr 2 ) + N c (51 + 44vr 2 - 216C(3)) + 9C F (3 - 4tt 2 + 48C(3)) 

I — 



K+(x) 
K+(x) 
Kt(x) 



2C F (-2C F + N c )x 
1 + x 

C F {9C F -llN c + 2n f )x 
3(x- 1) 

C f Nqx 



1 — X 

F" 



AC%x 



1 — x 
1 



K 6 + (x) = ~C F (10n f + N c (-67 + 3tt 2 )) 

y 

K i ( x ) = lc F (Wn f + N c (-Q7 + 3tt 2 )), 



(30) 



and for the (— ) case 

Kr(x) = ±C F (-2n f (3 + 4** 

_ 2C F (+2C F - N c )x 
2W " 1 + x 



N c (51 + Mir 2 - 216C(3)) + yC F (3 - 4vr 2 + 48C(3)) 
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KZ(x) 



C F (9C F - HJVo + 2n f )x 
3(x- 1) 

C F N c x 
1 — x 



T r / \ A:C rpX 

= -ic F (10 n/ + iV c (-67 + 3ir 2 )) 
y 

Kj(x) = --0(10/1/ - 18C F (x - 1) + iV c (-76 + 3tt 2 + 9x)). 



(31) 



The terms containing similar distribution (such as "+" distributions and 5 functions) have 
been combined together in order to speed-up the computation of the recursion relations. 



5 Comparisons among Moments 

It is particularly instructing to illustrate here briefly the relation between the Mellin 
moments of the parton distributions, which evolve with standard ordinary differential 
equations, and those of the arbitrary coefficient A n (x) and B n (x) which characterize 
Rossi's expansion up to next-to leading order (NLO). This relation, as we are going to 
show, involves a resummation of the ordinary moments of the parton distributions. 

Specifically, here we will be dealing with the relation between the Mellin moments of 
the coefficients appearing in the expansion 

A(N) = f 1 dxx N ~ 1 A(x) 



B(N) = f 1 dxx N ' 1 B(x) 
Jo 

(32) 

and those of the distributions 

A T q^(N,Q 2 ) = f 1 dxx N ~ l A T q { - ± \x,Q 2 )). (33) 
J o 

For this purpose we recall that the general (non-singlet) solution to NLO for the latter 
moments is given by 

(n (c,2\\-^tP^{N)I^ 

A T q ± (N,Q 2 ) = K{QlQ\N)[^J-) A T q±(N, Q 2 ) 

with the input distributions A T q±(Ql) at the input scale Q . We also have set 



(34) 



In the expressions above we have introduced the corresponding moments for the LO and 
NLO kernels (A T P$' N , A T P^f). 

The relation between the moments of the coefficients of the non-singlet x— space ex- 
pansion and those of the parton distributions at any Q, as expressed by eq. (JHIJl can be 
easily written down 

A n (N) + a s B n (N) = A T q ± (N, Q 2 )K(Q , Q, N) f ~ 2 ^ wW ) • (35) 
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As a check of this expression, notice that the initial condition is easily obtained from 
(|3ol) setting Q — > Q , n — > 0, thereby obtaining 



A^(N) +a s B^(N) = A T q ± (N,Qi), (36) 

which can be solved with A$ S (N) = A T q±(N,Q 2 ) and B^ S (N) = 0. 

It is then evident that the expansion (f20|) involves a resummation of the logarithmic 
contributions, as shown in eq. (pTH|) . 

In the singlet sector we can work out a similar relation both to LO 



with 



ei = -^-—(P i0) (N)-\ 2 l 



Aj — A 2 
1 



e 2 



-p(°)(N)+X 1 l 



A 2 — X\ 

Ai,2 = \ (p£> (N) + (iV) ± y (pj 9 0) (N) - pj? (iV)) 2 + 4P ( £ ) (iV)Pj? (iV)) , (38) 



2 

and to NLO 
where 



A„(A0 + a s P n (A0 = Xl (^^ +X 2 (-^) , (39) 



a f-2 e 2 Rei 
Xi = ei + — — eiRei + 



2vr V/5o A 1 -A 2 -/3 /2 / 

a f-2 eiRe 2 
X2 = e 2 + — — e 2 Re 2 + 



with 



2vr V/?o A 2 -A 1 -/3 /2 / 



R = P« (iV) - ApW (jv) . (41) 
2/?o 

We remark, if not obvious, that A n (N) and B n (N), pW(N), P W (N), in this case, are all 
2-by-2 singlet matrices. 



(40) 



6 Initial conditions 

As input distributions in the unpolarized case, we have used the models of Ref.pT], valid 
to NLO in the MS scheme at a scale Q 2 = 0.40 GeV 2 

x{u-u){x,Ql) = 0.632a; a43 (l -x) 3 ' 09 (l + 18.2x) 

x(d-d)(x,Q 2 ) = 0.624(1 - xf-°x{u -u){x,Ql) 

x(d-u)(x,Ql) = 0.20a; a43 (l -x) 12 - 4 (l - 13.3v^ + 60.0x) 

x(u + d)(x,Q 2 ) = 1.24x a20 (l -x) 8 - 5 (l -2.3v^ + 5.7x) 

xg(x,Ql) = 20.80x L6 (l-x) 4 - 1 (42) 

and xqi(x, Q 2 ) = xql(x, Qfc) = for % = s, c, b, t. 
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Following [8J, we have related the unpolarized input distribution to the longitudinally 
polarized ones 



xAu(x,Ql) = 1.019x a52 (l -xf 12 xu(x,Q 2 ) 



01 



xAd(x,Q 2 ) = -0.669x a4d xrf(x,Q 

xAu(x,Q 2 ) = -0.272x°- 38 xu{x,Ql) 

xAd(x, Ql) = xAu(x, Qq) 

xAg(x,Ql) = lA19x 1A3 (l-x)°- 15 xg(x,Ql) (43) 

and xAqi(x,Ql) = xAqi(x,Qo) = for q,i = s,c,b,t. Being the transversity distribution 
experimentally unknown, following [9], we assume the saturation of Soffer's inequality 



2^ 



xA Tqi (x, Ql) = Mx,Ql)+xA qi (x,Q 0J _ (44) 

7 Names of the input parameters, variables and of the 
output files 

7.1 Notations 

gluons, g g 

1-6 quarks, % sorted by their mass valuesfw, d, s, c, b, t) u.d.s.c.b.t 

7-12 antiquarks, qj au,ad,as,ac,ab,at 

13-18 q\ ' um , dm , sm , cm , bm , tm 

19-24 Xi (unpolarized and longitudinally polarized cases) Cu,Cd,Cs , Cc , Cb ,Ct 

qf^ (transversely polarized case) CUjCd^s ,Cc ,Cb,Ct 

25 qr(+> qp 



7.2 Input parameters and variables 

process unpolarized 

1 longitudinally polarized 

2 transversely polarized 
spacing 1 linear 

2 logarithmic 

GRID_PTS Number of points in the grid 

MGP Number of Gaussian points, ric 

ITERATIONS Number of terms in the sum (J2DJ) 

extension Extension of the output files 
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betaO 
betal 
alphal 
alpha2 



step 

lstep 

X[i] 

XG[i] [j] 
WG[i] [j] 
nf, Nf 
n_evol 



lambda [i] 
A[i] [j] [k] 
B[i] [j] [k] 



Q[i] 




7.3 Output files 



The generic name of an output file is: Xni.ext, where 

X is U in the unpolarized case, L in the longitudinally polarized case and T in the trans- 
versely polarized case; 

n is a progressive number that indicates the scale Q 2 at a given stage: n = refers to 
the initial scale, the highest value of n refers to the final scale and the intermediate 
values refer to the quarks production thresholds (1 for charm, 2 for bottom and 3 
for top); 

i is the identifier of the distribution, reported in the third column of the table in subsec- 
tion ED 

ext is an extension chosen by the user. 

8 Description of the program 
8.1 Main program 

At run time, the program asks the user to select a linear or a logarithmic spacing for the 
x-axis. The logarithmic spacing is useful in order to analyze the small-x behavior. Then 
the program stores as external variables the grid points X{ and, for each of them, calls the 
function gauleg which computes the Gaussian points and weights Wij corresponding 
to the integration range [xi, 1], with < j < tiq — 1. After that, the user is asked to enter 
the type of process, the final value of Q and an extension for the names of the output 
files. At this point the program computes the initial values of the parton distributions 
for gluons, up, down, anti-up and anti-down (see Section EJ) at the grid points and stores 
them in the arrays A[i] [0] [k] (see (f2"lj) ). setting to zero the initial distributions of the 
heavier quarks. 

The evolution is done in the various regions of the evolutions, all characterized by 
a specific flavour number. Each new flavour comes into play only when the scale Q 
reaches the corresponding quark mass. In that case rif is increased by 1 everywhere in 
the program. The recurrence relations (|22|) and (J2*3*jl are then solved iteratively for both 
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the nonsinglet and the singlet sector, and at the end of each energy step the evolved 
distributions are reconstructed via the relation (ffill) . The distributions computed in this 
way become the initial conditions for the subsequent step. The numerical values of the 
distributions at the end of each energy step are printed to files. 

8.2 Function writef ile 

void writef ile (double *A,char *f ilename) ; 

This function creates a file, whose name is contained in the string *f ilename, with 
an output characterized by two columns of data: the left column contains all the values 
of the grid points Xi and the right one the corresponding values of the array A [i] . 

8.3 Function alpha_s 

double alpha_s (double Q, double betaO, double betal, double lambda); 

Given the energy scale Q, the first two terms of the perturbative expansion of the 
/3-function betaO and betal and the value lambda of alpha_s returns the two-loop 

running of the coupling constant, using the formula . 

8.4 Function gauleg 

void gauleg(double xl,double x2,double x[],double w[],int n) ; 

This function is taken from [10] with just some minor changes. Given the lower and 
upper limits of integration xl and x2, and given n, gauleg returns arrays x[0, . . . 5 n-l] 
and w [0 , . . . ,n-l] of length n, containing the abscissas and weights of the Gauss-Legendre 
n-point quadrature formula. 

8.5 Function interp 

double interp (double *A, double x) ; 

Given an array A, representing a function known only at the grid points, and a number 
x in the interval [0, 1], interp returns the linear interpolation of the function at the point 
x. 

8.6 Function IntGL 

double IntGL(int i,double kernel(double z), double *A) ; 

Given an integer i (corresponding to a grid point x^), a one variable function kernel (z) 
and an array A, representing a function g(x) known at the grid points, IntGL returns the 
result of the integral 




(45) 



computed by the Gauss-Legendre technique. 



8.7 Function IntPD 



double IntGL (int i, double *A) ; 
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Given an integer i, to which it corresponds a grid point X{, and an array A, representing 
a function f(x) known at the grid points, IntGL returns the result of the convolution 



® /(»<)=/ + f{Xi)log{l - Xi), (46) 



[i — Xi) + Jxi y y — Xi 

computed by the Gauss-Legendre technique. 



8.8 Function S2 

double S2 (double z) ; 

This function evaluates the Spence function 5*2(2) using the expansion 

1 tt 2 00 (— 7 \ n 

S 2 (z) = logzlog(l -z)-- (\ogzf + - + E Mr" ( 47 ) 

4 iz n =l n 

arrested at the 50th order. 



8.9 Function fact 

double fact(int n) ; 

This function returns the factorial n! 



8.10 Initial distributions 

double xuv (double x) ; 

double xdv (double x) ; 

double xdbmub (double x) ; 

double xubpdb (double x) ; 

double xg (double x) ; 

double xu (double x) ; 

double xubar (double x) ; 

double xd (double x) ; 

double xdbar (double x) ; 

double xDg (double x) ; 

double xDu (double x) ; 

double xDubar (double x) ; 

double xDd (double x) ; 

double xDdbar (double x) ; 

Given the Bjorken variable x, these functions return the initial distributions at the 
input scale (see Section EJ). 



8.11 Regular part of the kernels 

double PONS (double z) ; 

double P0qq(double z) ; 

double POqg (double z) ; 

double P0gq(double z) ; 

double POgg (double z) ; 

double PlNSm(double z) ; 

double PINSp (double z) ; 

double Plqq(double z) ; 
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double 


Plqg(.double z 


) , 


double 


Plgq(.double z 


) , 


double 


Plgg (double z 


) ; 


double 


DPONS (double 


z) ; 


double 


DPOqq (.double 


z; ; 


double 


DPOqg (double 


z) ; 


double 


DPOgq (.double 


z; ; 


double 


DPOgg (double 


z) ; 


double 


DP INSm (double 


z) 


double 


TN T>1 A TIT/") / 1 IT 

DPIMSp (double 


z) 


double 


DP lqq (.double 


z; ; 


double 


DPlqg (double 


z) ; 


double 


DPlgq(double 


z) ; 


double 


DPlgg (double 


z) ; 


double 


tPO (double z) 


j 


double 


tPlm(double z 


) ; 


double 


tPlp (double z 


) ; 



Given the Bjorken variable z, these functions return the part of the Altarelli-Parisi 
kernels that does not contain singularities. 

9 Running the code 

In the plots shown in this paper we have divided the interval [0, 1] of the Bjorken variable 
x in 500 subintervals (GRID_PTS=501), 30 Gaussian points (NGP=1), and we have retained 

10 terms in the sum J2H) (ITERATI0NS=10). In the figures and d the flag spacing has 
been set to 2, in order to have a logarithmically spaced grid. This feature turns useful 
if one intends to analyze the small-x behavior. We have tested our implementation in a 
detailed study of Soffer's inequality up to NLO [T^] . 

0.9 
0.8 
0.7 
0.6 

s 05 

x 0.4 
0.3 
0.2 
0.1 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x 

Figure 1: Evolution of the unpolarized quark up distribution xu versus x at various Q 
values. 



10 Conclusions 

We have illustrated and documented a method for solving in a rather fast way the NLO 
evolution equations for the parton distributions at leading twist. The advantages of the 
method compared to other implementations based on the inversion of the Mellin moments, 



Q=0.632 GeV 




15 




as usually done in the case of QCD, are rather evident. We have also shown how Rossi's 
ansatz, originally formulated in the case of the photon structure function, relates to the 
solution of DGLAP equations formulated in terms of moments. The running time of the 
implementation is truly modest, even for a large number of iterations, and allows to get a 
very good accuracy. We hope to return with a discussion of the supersymmetric extension 
of the method and other applications in the near future. 
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A Unpolarized kernels 

The reference for the unpolarized kernels is |2j, rearranged for our purposes. We remind 
that the plus distribution is defined by 
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Figure 5: Evolution of the unpolarized gluon distribution xg versus x at various Q values. A 



and the Spence function is 



1 7T 2 
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Zi D 



where the dilogarithm is defined by 
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Figure 6: Evolution of the longitudinally polarized quark up distribution xAu versus x 
at various Q values. 
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Figure 7: As in figure (0), but now the x-axis is in logarithmic scale, to show the small-x 
behavior. 
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Figure 8: Evolution of xAd versus x at various Q values. 
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Figure 9: Evolution of xAs = xAs versus x at various Q values. 
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Figure 10: Evolution of the longitudinally polarized antiquark up distribution xAu versus 
x at various Q values. 
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Figure 11: Evolution of the longitudinally polarized gluon distribution xAg versus x at 
various Q values. 
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Figure 12: As in figure (jUJ), but now the x-axis is in logarithmic scale. 
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Figure 13: Evolution of the transversely polarized quark up distribution xA F u versus x 
at various Q values. 
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Figure 14: Evolution of xAxd versus x at various Q values. 
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Figure 15: Evolution of the transversely polarized antiquark up distribution xAtU 
x at various Q values. 
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B Longitudinally polarized kernels 

The reference for the longitudinally polarized kernels is [TT| . 
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+ 



+ 



+ 



' 4C F 7)x(x 2 - 1) + N%(1 + x(6 + x(2 + (x - 8)x))) 

2(1 -x)x 
'2iV 2 (l -x(2-2x-x 3 )) 



log X 



(x — l)x 
N^il + x + 3x 2 + x 3 ) 

X 



I logxlog(l — x) 

S 2 (x) 



26 



+ 



N, 



N C (Q7 - 3tt 2 ) - 207) j 
+ {AS(| + 3C(3))-^T / -^>} ff (l-x) 



(1-*) 
47Vc7) ' 



P. 



(i) 



(x) 



18a; 



[6Cf7)(:e - l)(x(7 + lOx) - 2) 

+iv£(2(70 - 3x)x - 3tt 2 (1 - a; + 3a; 2 - a; 3 )) 
+2N c T f (x(37 + x(23x - 57)) - 23)]} 

+ <2C F T f (l - 3ar) - iV 2 f 9 - 13x + ^f-) | log a; 



+ {^(l-x + 3x 2 



+ 



2x 

x 



X 



log X 



(a; 3 — 3x 2 + x — | log a; log(l 



— x 



r_iVg(l + »(2 + 2, + *»))} 
1 x(l+x) J 



D Transversely polarized kernels 

The reference for the transversely polarized kernels is |12j . 



A T P 



(i) 



9 



207) - 18C F (x - 1) + N c (9x - 76 + 3vr 2 



+ 



+ 
+ 



'C F (9C F -llN c + 4T f )x' 
3(1-1) . 

Cf^OSI 2 

— f log X 

1 — X J 

4C 2 x 



log a; 



' 1 logxlog(l — x) 



x — 1 

' 2C F (2C F - N c )x 



+ 



9 

cv 

72 



l + x 

N c (3n 2 - 67) + 2QT 



(1-x). 



iV c (51 + 44tt 2 - 216C(3)) - 47) (3 + 4vr 2 
+9C F (3 - 4vr 2 + 48C(3))] } 5(1 - x) 



A P (i) ... 1°^ 
a tP ns + - < — 



9 



+ 



iV);(37r 2 - 67) + 207) 

f CV(9CV~ lliV c + 47)) 
\ 3(x- 1) 



log a; 



27 



■ C F N c x \ 2 
+ i — f log x 



+ 



1 — X 
AC 2 F x \ 
x-lj 

2C F (N C - 2C F )x 



log x log (1 — x) 

S 2 (x) 



9 

C F 



1 + X 

N c (3n 2 - 67) + 2QT- 



1 1 <l-x) + 

2 11 (?/■/ Q\\ /IT to I A~3\ 



N c (51 + 44tt 2 - 216C(3)) - 4T/(3 + An 
+9C F {3 - 4tt 2 + 48C(3))] } 5(1 - x) 



(85) 
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